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ABSTRACT 

We introduce the Phoenix Project, a set of ACDM simulations of the dark matter component 

■ of nine rich galaxy clusters. Each cluster is simulated at least at two different numerical reso- 
Qh lutions. For eight of them, the highest resolution corresponds to ~ 1 30 million particles within 

the virial radius, while for one this number is over one billion. We study the structure and sub- 
£^ ■ structure of these systems and contrast them with six galaxy-sized dark matter haloes from 

^ | the Aquarius Project, simulated at comparable resolution. This comparison highlights the ap- 

proximate mass invariance of CDM halo structure and substructure. We find little difference 
in the spherically-averaged mass, pseudo-phase-space density, and velocity anisotropy pro- 
files of Aquarius and Phoenix haloes. When scaled to the virial properties of the host halo, 
the abundance and radial distribution of subhaloes are also very similar, despite the fact that 
1 Aquarius and Phoenix haloes differ by roughly three decades in virial mass. The most notable 

t^j- , difference is that cluster haloes have been assembled more recently and are thus significantly 

■ less relaxed than galaxy haloes, which leads to decreased regularity, increased halo-to-halo 

scatter and sizable deviations from the mean trends. This accentuates the effects of the strong 
asphericity of individual clusters on surface density profiles, which may vary by up to a factor 

■ of three at a given radius, depending on projection. The high apparent concentration reported 
I for some strong-lensing clusters might very well reflect these effects. A more recent assem- 
bly also explains why substructure in some Phoenix haloes is slightly more abundant than 
in Aquarius, especially in the inner regions. Resolved subhaloes nevertheless contribute only 
1 1 ±3% of the virial mass in Phoenix clusters. Together, the Phoenix and Aquarius simulation 
series provide a detailed and comprehensive prediction of the cold dark matter distribution in 

5^ 1 galaxies and clusters when the effects of baryons can be neglected. 



Key words: methods: N-body simulations - methods: numerical -dark matter - galaxies: 
haloes - galaxiesxlustering 



1 INTRODUCTION 

The past two decades have witnessed the emergence of a paradigm 
for the origin of structure in the Universe. There is now strong ev- 
idence that the dominant forms of the matter-energy content are a 
combination of a mysterious form of "dark energy" that governs 
the late expansion of the Universe, and "dark matter" made up of 
some kind of non-baryonic, weakly interacting elementary particle 
left over from the Big Bang. Although the exact nature of the dark 
matter particle is unknown, astrophysical clues to its identity may 
be gained by studying its clustering properties on different scales. 
Considerable effort has been devoted to this task, and has led to the 
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crafting of detailed theoretical predictions, especially for the case 
of particles with negligible thermal velocity, the cornerstone of the 
popular "cold dark matter" (CDM) theory. As a result, we now un- 
derstand fairly well: (i) the statistics of CDM clustering on large 
scales and its dependence on the cosm ological parameters (e.g., 
I Jenkins et all 1 19981 ; ISpringel et al]|2006l) ; (ii) the dyn amics of its 
incorp oration into non-linear units ("haloes") (see, e.g. JWang et al.l 
1201 ll and references therein); and, at least empirically, (iii) its spa- 
tial d i stribution within such virialized struct ures (e.g.. lFrenk et al.l 
1 19851 ; INavarro et al Jl 1 996L [l997l . 120041 

Progress in this field has been guided by N-body sim- 
ulations of ever increasing numerical resolution a nd dynamic 
range (e.g. lFrenketal.1 19851; iNavarro et al .11 19971: iMoore et al.l 



1 19991; Ijing & Sutol 12003; INavarro et al.l |2004| ; iGao et alj 12004b: 



Gao et al. 



Diemand et alj2004<i 20071: Gao et al 


20081: ISpringel et al. 2008a; 


Stadel et al. 20091; Navarro et alj 201C 


). These simulations are es- 



sential to investigate highly non-linear scales such as the haloes 
of individual galaxies and galaxy groupings, where simple analyt- 
ical approximations fail. A few key properties of CDM haloes are 
now widely agreed upon, at least when the effects of baryons are 
neglected: (a) the presence of a central density "cusp"; (b) strong 
deviations from spherical symmetry; (c) a remarkable similarity in 
the shape of the mass profiles; and (d) the presence of abundant 
substructure in the form of self-bound "subhaloes". 

On the scale of individual galaxies, these key predictions have 
been confirmed and refined by the lates t simulation work, i n par- 
ticular the Via Lactea simulation series (Diemand et al. 200^), the 



GHALO simulation (Stadel eta 1. 2009) and the Aquarius Project 
( Sprin gel et alj2008bl lal: lNavarro et alj 2010). For example, the cen- 



tral density cusp is now accepted to be shallower than hypothesized 
in some earlier work and mass profiles have been shown to be only 
approximately self-similar. Further, it is now clear that although 
subhaloes are subdominant in terms of total mass, they are still 
dense and abundant enough to dominate the dark matter annihi- 
lation radiation from a halo. 

As shown bv lSpringel eTail d2008ah . the latter statement re- 
quires a detailed characterization of the substructure, including the 
internal properties of the subhaloes, their mass function, and their 
spatial distribution within the main halo. The Aquarius Project has 
provided compelling, if mainly empirical, guidance to each of these 
issues in the case of haloes similar to that of the Milky Way. For ex- 
ample, the subhalo mass function is well approximated by a power 
law, dN/dM <x M~ 19 , wi th a normalisation , in scaled units, weakly 
dependent on halo mass dGao et alj|207l^ . In addition, subhaloes 
tend to avoid the central region of the main halo and are more 
prevalent in the outer regions. Interestingly, their spatial distribu- 
tion appears independent of subhalo mass, a result that, if generally 
applicable, simplifies substantially the characterization of substruc- 
ture. Finally, the internal structure of subhaloes obeys scaling laws 
similar to those of haloes in isolation but slightly modified by the 
effects of the tidal field of the main halo: subhaloes are "denser", 
reaching their peak circular velocity at radii roughly half that of 
their isolated counterparts. 

Galaxy clusters are a promising venue for testing these pre- 
dictions. The central cusp, for example, can be constrained by 
combining measurements of the stellar kinematics of the central 
galaxy with a lensing a nalysis of radial a n d tangential "arcs" near 
the c l uster centre (e.g.. ISand et alj|2002l. [20041 ; iMeneghetti et al] 
120071 : iNewman et"aT]|2009t IZitrin et al] |201 lh . Outside the very 
centre, the clu ster mass profile can be measured through weak lens- 



ing ( see, e.g., lOkabe et al.ll20ld : bguri et aljboi ll : llJmetsu et al] 



2011), X-ray stu dies of the hot intracluster medium (ICM; e.g, 



Buote et alj 2007), and, more recently, through the ICM Sunyaev- 
ZeTdovich effect on the cosmic microwave background (see, e.g., 
iGralla etal]|201 J). In many cases, including s ubstructure seems 
required in order to obtain acceptable fits (e.g. iMao & Schneider! 
ll998l:lMao et al.l2004IXu et alj2009l : lNataraian et alj200ll2009n T 
implying that it should be possible to contrast observations directly 
with the CDM substructure predicted by simulations. 

Such endeavour has so far been hindered by the lack of ultra- 
high-resolution dark matter simulations of galaxy clusters compa- 
rable to the Aquarius series. Indeed, the highest-resolution galaxy 
cluster simulations published to date have at most of order a few 
million particles with i n the virial rad i us (e.g . Jing & Sutdl2000; 
ISpringel et all 1200 lj ; iDiemand et all l2004a[ ; iReedet al] l2005h . 
roughly one thousand times fewer than the best resolved Aquar- 



ius halo. None of these cluster simulations are thus able to address 
conclusively issues such as the structure of the central cusp or the 
properties of cluster substructure. 

Although it may be tempting to appeal to the nearly self- 
similar nature of CDM haloes to extrapolate the Aquarius results 
to cluster scales, it is unclear what systematic uncertainties might 
be introduced through such extrapolation. Clusters are rare, dy- 
namically young objects up to one thousand times more massive 
than individual galaxies. They thus trace scales where the CDM 
power spectrum differs qualitatively from that of galaxies. Preci- 
sion work demands that the near self-similarity of dark haloes be 
scrutinized directly in order to provide definitive predictions for the 
CDM paradigm on cluster scales. 

To this aim, we have carried out a suite of simulations de- 
signed to address these issues in detail. The Phoenix Project fol- 
lows the design of the Aquarius Project and consists of zoomed-in 
resimulations of individual galaxy clusters drawn from a cosmo- 
logically representative volume. The simulations follow only the 
dark matter component of each cluster, and include the first simu- 
lation of a cluster-sized halo with more than one billion particles 
within the virial radius. Like the Aquarius Project on galaxy scales, 
the large dynamic range of these simulations allows us to probe 
not only the innermost regions of cluster haloes and thus the struc- 
ture of the central cusp, but also the statistics, internal structure, 
and spatial distribution of cluster substructure over a mass range 
spanning seven decades. 

Our paper is organized as follows. In Section [2] we describe 
our numerical techniques and introduce the simulation set. In Sec. [3] 
and Sec.|4]we discuss, respectively, the density profile and substruc- 
ture properties of Phoenix haloes and compare them with those of 
Aquarius. Sec.[5]summarizes our main conclusions. 



2 THE SIMULATIONS 

The Phoenix Project consists of a series of simulations of 9 differ- 
ent galaxy clusters with masses exceeding 5 x 10 14 hT x Mq. These 
clusters were selected from a large cosmological volume and resim- 
ulated individually at varying resolution. Details of the resimulation 
procedure are given below. 



2.1 Cosmology 

All the simulations reported here ado pt the cosmologica l param- 
eters of the Millennium Simulation l lSpringel et alJuOOa) : O.^ = 
0.25, flyv = 0.75, (78 = 0.9, n s = 1, and a present-day value of the 
Hubble constant H = 100/ikm Mpc~' = 73 km s -1 Mpc~' . 
This is also the s et of cosmological p arameters adopted for the 
Aquarius project ( ISpringel et alj|2008bl) . which targeted haloes a 
thousand times less massive. Although th ey are inconsistent with 
the latest CMB data (Kom atsu et al.ll20l"lh the differences are not 
large (the main difference is that a lower value of =0.81 is now 
preferred) and they are expected to affect only the abundance of 
cluster hal oes rather than the ir detailed structure and substructure 
properties dWang et al.ll2012h . This choice also has the advantage 
that any difference between Aquarius and Phoenix haloes can be 
traced to the different mass scales and not to variations in the cos- 
mological model. 
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Name 




M 2 oo 


^200 


N200 


E 


F conv 






Iff 1V10J 


III 1V1UC 






r/7 _1 knrl 

It JvUL. 


Ph-A-1 


6.355 x 10 5 


6.560 x 10 14 


1.413 


1,032,269.120 


0.15 


1.2 


Ph-A-2 


5 084 y 1 D 6 


6 570 x 10 14 
u. j 1 \j a iu 


1.414 


1 7Q 935 47? 


32 


2 7 


Ph-A-3 


1.716 x 10 7 


6.566 x 10 14 


1.413 


38,261,560 


0.7 


4.2 


Ph-A-4 


1.373 x 10 s 


6.593 x 10 14 


1.415 


4,802,516 


2.8 


9.4 


Ph-B-2 


6.127 x 10 s 


8.255 x 10 14 


1.526 


134,718,112 


0.32 


3.0 


Ph-B-4 


1.656 x 10 s 


8.209 x 10 14 


1.522 


4,956,688 


2.8 


10.7 


Ph-C-2 


4.605 x 10 s 


5.495 x 10 14 


1.386 


119,324,008 


0.32 


2.6 


Ph-C-4 


1.182 x 10 8 


5.549 x 10 14 


1.383 


4,696,046 


2.8 


9.2 


Ph-D-2 


4.721 x 10 s 


6.191 x 10 14 


1.386 


130,529,200 


0.32 


2.7 


Ph-D-4 


1.373 x 10 s 


6.162 x 10 14 


1.384 


4,488,330 


2.8 


9.4 


Ph-E-2 


4.425 x 10 s 


5.969 x 10 14 


1.369 


130,529,200 


0.32 


2.4 


Ph-E-4 


1.017 x 10 8 


5.923 x 10 14 


1.366 


5,824,375 


2.8 


8.4 


Ph-F-2 


4.425 x 10 s 


7.997 x 10 14 


1.509 


129,221,216 


0.32 


2.8 


Ph-F-4 


1.682 x 10 8 


8.039 x 10 14 


1.512 


4,779,008 


2.8 


10.3 


Ph-G-2 


8.599 x 10 s 


1.150 x 10' 5 


1.704 


133,730,958 


0.32 


3.2 


Ph-G-4 


2.907 x 1() 8 


1.148 x 10 15 


1.703 


3,949,310 


2.8 


13.1 


Ph-H-2 


8.600 x 10 s 


1.136 x 10' 5 


1.686 


129,488,456 


0.32 


2.9 


Ph-H-4 


2.502 x 10 8 


1.150 x 10 15 


1.686 


4,456,720 


2.8 


11.8 


Ph-I-2 


1.841 x 10 7 


2.411 x 10' 5 


2.185 


131,845,620 


0.32 


2.9 


Ph-I-4 


4.559 x 10 8 


2.427 x 10 15 


2.181 


5,289,259 


2.8 


14.2 



Table 1. Basic parameters of the Phoenix simulations. Each of the nine haloes is labelled as Ph-X-N, where the letter X (from A to I) identifies each halo, and 
N, which runs from 1 to 4, refers to the numerical resolution (1 is highest). The parameter m v gives the particle mass in the high-resolution region that includes 
the cluster; M200 is the virial mass of the halo; noo is the corresponding virial radius; and /V200 states the number of particles inside (-200- The parameter £ is 
the Plummer-equivalent gravitational softening length, so that pairwise interactions are fully Newtonian when separated by a distance greater than 2.8 E. The 
last column lists the "convergence radius", r conv , outside which the circular velocity is expected to converge to better than 10%. 



2.2 Cluster Sample and Resimulations 

The Phoenix cluster sample is selected for resimulation from the 
Millennium Simulation friends-of-friends group catalog at z = 0. 
Six clusters were selected at random from the 72 systems with 
viriari mass in the range 5 < M 2 oo/1O 14 A _1 M < 10. In order 
to sample the tail of rare rich clusters three further Phoenix clus- 
ters were selected from the nine Millennium halos which have 
M 2 oo > W l5 hr l M®. 

The initial conditions for resimulation were set up using a 
procedure analogo us to that used for the Aquarius haloes and de - 
scribed in detail bv lPower etal] J2003I) and lSpringel et al.l d2008al) . 
The only difference is that the initial displacements and velocities 
were computed using second-or der Lagrangian perturbation theory, 
as described by I Jenkins! JioToh . All nine haloes were resimulated 
at least twice using different numerical resolution (level 2 and level 
4, respectively). At level 2 each cluster has between 120 and 135 
million particles within the virial radius; at level 4 each system is 
made up of 4 to 6 million particles. 

We have selected one of the clusters (Ph-A) for a numerical 
resolution study and have carried out an extra level-3 run (with 
roughly 40 million particles within ^200) an d a flagship level- 1 
run, where we followed 4.05 billion high-resolution particles in 



1 We define the virial radius of a cluster, rim, as that of a sphere of mean 
density 200 times the critical density for closure; p cr i t = 3Hq /8jiG. The 
virial radius defines implicitly the virial mass of a cluster, M200. and its 
virial velocity, V200 = \/ 'GM 2 oo / >200 ■ 



total, 1.03 billion of which are found within ^200 at z = 0. For 
ease of reference we label the runs using the convention Ph-X- 
N, where X is a letter from A to I that identifies each individual 
cluster and N is a number from 1 to 4 that specifies the resolution 
level. The simulation parameters are summarized in Table [T] We 
ha ve used for all runs the P-Gadget-3 code, a version of Gadget- 
2 ^Springel et al.ll2005T ) especially optimized for zoomed-in cos- 
mological resimulations in distributed-memory massively-parallel 
comput ers. The code is iden tical to that used for the Aquarius 
Project (Sprin gel etal] [2008b). The simulations were carried out 
on Deepcomp 7000 at the Supercomputer Center of the Chinese 
Academy of Science. The largest simulation, Ph-A-1, used 3 Tbs 
of memory on 1024 cores and took about 1.9 million CPU hours. 
The initial conditions were generated at the Institute for Computa- 
tional Cosmology (Durham University). 

The gravitational soften ing of each run wa s chosen following 
the "optimal" prescription of Power et al.l d2003h . It is kept fixed in 
comoving coordinates throughout each run and is listed in Table 
Q] Our highest-resolution run (Ph-A-1) has a nominal (Plummer- 
equivalent) spatial resolution of just 150/i~' pc. 

Haloes are identified in each run using the friends-of-friends 
(FOF) group finding algorithm with linking l ength set to 20% of 
the mean interparticle separation l lDavis etalll 19851). Substructur e 
within FOF haloes is identified by SUBFIND ( Spri ngel et al.l2001bl) . 
a groupfinding algorithm that searches recursively for self-bound 
subhaloes. Both FOF and SUBFIND have been integrated within P- 
Gadget-3 and are run on-the-fly each time a simulation snapshot is 
created. 
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Figure 1. Images of cluster Ph-A at four different numerical resolutions. Each panel projects a cubic volume 5h Mpc on a side. The brightness of each 
image pixel is proportional to the logarithm of the square of the dark matter density projected along the line of sight, and the hue encodes the local velocity 
dispersion density-weighted along the line of sight (see text for details). This rendering choice highlights the presence of substructure which, although abundant, 
contributes less than about 10% of the total mass within the virial radius. 



We have stored for each run 72 snapshots uniformly spaced in 
log 1() a, starting at a = 0.017 (a = 1/(1 +z) is the expansion factor). 
The initial conditions are set at Zi n i t = 63 for our level 4 and at 
Zinit = 79 for the rest. The large number of outputs is designed 
to allow us in future work to implement semi-analytic models of 
galaxy formation in order to follow the evolution of the baryonic 
component of galaxies within rich clusters. 

We list the basic structural parameters of Phoenix clusters at 
redshift z = in Table [2] These include the peak circular velocity, 
Vmax. an d the radius, r max , at which it is reached; the half-mass for- 
mation redshift, Zh, when the main progenitor first reaches half the 
final halo mass; the conc entration parameters , c and eg, obtained 
from the best-fit NFW jNavarro et aHll996t 1 19971) and Einasto 
(1965) profiles, respectively; the figure of merit, <2nfw and Qe, as- 



sociated with each of those fits; and the Einasto "shape" parameter 
a. (See the Appendix for definitions corresponding to these fitting 
formulae and for details on the profile-fitting procedure.) A' su b is 
the total number of subhaloes with more than 20 particles identi- 
fied by SUB FIND inside r2oo; /sub ls the total mass contributed by 
these subhaloes, expressed as a fraction of the virial mass. 



3 THE STRUCTURE OF PHOENIX CLUSTERS 

We shall focus our analysis on the properties of Phoenix clusters 
at z = 0. Figure[T] sho ws Ph-A at the four different numerical reso- 
lutions. As in Sprin geletal"1 (l2008a). this and other cluster images 
are constructed so that the brightness of each pixel is proportional 
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Figure 2. The inner ~ l/i 1 Mpc of Ph-A-1. Color coding is as in Fig. [T] This figure illustrates clearly the strong asphericity of the halo; the presence of 
several nested levels of substructure, and the tendency of subhaloes to avoid the halo centre. 



to the logarithm of the square of the dark matter density projected 
along the line of sight, 

S{x,y) = ftf K {r)di (1) 

while the color hue encodes the mean dark matter velocity disper- 
sion, 

o{x,y) = -^—r J o loc (r) pf oc (r) dz (2) 

Here the local dark matter density, Pi oc (r), and the local velocity 
dispersion, Oi oc (r), are estimated using an SPH kernel interpolation 
scheme. 

Figure [7] shows that the main result of increasing the number 
of particles is the ability to resolve larger numbers of subhaloes. On 
the other hand, the main properties of the cluster, such as its shape 



and orientation, the overall mass profile, and even the location of 
the largest subclumps remain invariant in all four Ph-A realizations. 

Fig. [2] is analogous to Fig. Q] but for the inner ~ IhT 1 Mpc 
of Ph-A-1 (our highest resolution run). This image highlights the 
strong asphericity of the halo, as well as the presence of several 
nested levels of substructure (i.e., subhaloes within subhaloes). It 
also shows that subhaloes tend to avoid the central regions. These 
charac teristics are shared with galaxy-sized haloes dSpringel et al.l 
2008a), and appear to be typical of CDM haloes on all mass scales. 

Fig.[3]is analogous to Fig.Q]but for all level-2 Phoenix haloes 
at z = 0. This figure shows that the main characteristics of Ph-A de- 
scribed above are common to all Phoenix clusters: strong aspheric- 
ity; abundant substructure; and a marked difference between the 
spatial distribution of mass (which is highly concentrated) and that 
of subhaloes (which tend to avoid the central regions). 



6 Gao et al. 





Figure 3. As Fig.[TJ but for all level-2 Phoenix clusters at z = 0. Boxes are all 5h~ [ Mpc on a side. Note that the appearance of several Phoenix clusters is 
suggestive of a transient evolutionary stage, characterized by the presence of a number of undissolved substructure groupings. Ph-G-2 is a particularly good 
example of this irregular structure which may be traced to its recent assembly time; this cluster has acquired half its mass since z = 0. 18. 



Fig. [3] also highlights an important characteristic of cluster- 
sized dark matter haloes: the presence of "multiple centres" traced 
by groups of subhaloes, as well as the overall impression that many 
systems are in a transient, unrelaxed stage of their evolution. This 
is expected, given the late assembly of the clusters: Ph-G-2, for 
example, assembled half its final mass after z = 0.18; the median 
half-mass assembly redshift for all Phoenix clusters is just z = 0.56. 
Ph-A, on the other hand, appears relaxed; this cluster has the high- 
est formation redshift of our sample, Zh ~ 1.2. 

The late assembly and concomitant departures from equilib- 
rium are characteristics that set clusters apart from galaxy-sized 
haloes; for comparison, the median half-mass formation redshift of 
Aquarius haloes is z ~ 2. Table [2] lists two quantitative measures 
of departures from equilibrium: the fraction of mass in substruc- 
tures, / su b, and the offset, d D ff, between the centre of mass of the 



halo and the location of the potential minimum expressed in units 
of the virial radi us (for further discussion of these parameters see 
iNeto et al] [2007"). These correlate well with the formation redshift, 
Zh, and are significantly larger, on average, than in the galaxy-sized 
Aquarius haloes (see Tabled- 

3.1 Mass Profiles 

We explore in this section the spherically-averaged mass profiles of 
Phoenix haloes. We begin by using the four Ph-A realizations to as- 
sess the limitations introduced by finite numerical resolution. Fig. [4] 
shows the density profile, p(r), as well as the radial dependence of 
the logarithmic slo pe, y= — dlnp/dlnr , for Ph-A-1 throug h Ph-A - 
4. As discussed by iPower et"aTI d2003h and iNavarro et al.l J2010l) . 
the mass profiles of simulated haloes are robustly determined in 
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Figure 4. Left panel: Spherically-averaged density profile of halo Ph-A at z = 0. Different colors correspond to the four different resolution runs listed in 
Table [T] The panel on the left shows the density multiplied by r 2 in order t o enhance the dyn amic range of the plot. Each profile is shown with a thick 
line connecting filled circles from the "convergence radius", r conv , outwards (Power etal. 2003]). Thin curves extend the profiles inwards down to r = 2e, 
where 8 is the Plummer-equivalent gravitational softening length. Vertical dotted lines indicate, for each run, 2.8 e, the distance beyond which pairwise particle 
interactions are fully Newtonian. Note the excellent numerical convergence achieved for each simulation outside their r conv . An NFW profile with concentration 
c = 5.63 (thin dashed brown line) and an Einasto profile with a = 0.22 and eg = 5.59 (thin dashed magenta line) are also shown for comparison. Right panel: 
Logarithmic slope (y= — dlnp/dlnr) of the density profile as a function of radius. Colors and line types are the same as in the left panel. Note again the 
excellent convergence achieved in all runs at radii outside the convergence radius, r COIlv . 
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''max 


;h 




c 


Qe 


2nfw 


a 


Wsub 


/sub 


doff 




[fans" 1 ] 


[h^Mpc] 




















Ph-A-1 


1521.82 


0.55 


1.17 


5.59 


5.63 


0.037 


0.093 


0.215 


192,206 


0.080 


0.04 


Ph-A-2 


1527.24 


0.55 


1.17 


5.72 


5.96 


0.039 


0.075 


0.216 


26,896 


0.071 


0.04 


Ph-A-3 


1529.41 


0.56 


1.17 


5.69 


6.04 


0.038 


0.061 


0.218 


8,478 


0.062 


0.04 


Ph-A-4 


1538.88 


0.59 


1.17 


5.71 


6.14 


0.052 


0.063 


0.219 


1,049 


0.049 


0.04 


Ph-B-2 


1624.52 


0.53 


0.46 


4.41 


4.19 


0.127 


0.108 


0.235 


38,659 


0.108 


0.02 


Ph-B-4 


1623.12 


0.56 


0.46 


4.40 


4.06 


0.107 


0.117 


0.276 


1,657 


0.081 


0.02 


Ph-C-2 


1294.19 


0.65 


0.76 


4.27 


5.11 


0.077 


0.104 


0.181 


33,529 


0.114 


0.06 


Ph-C-4 


1310.19 


0.78 


0.76 


4.34 


4.72 


0.085 


0.112 


0.185 


1,489 


0.095 


0.06 


Ph-D-2 


1393.13 


0.68 


0.46 


3.88 


4.08 


0.122 


0.086 


0.205 


38,199 


0.124 


0.05 


Ph-D-4 


1436.10 


0.65 


0.46 


4.03 


4.34 


0.136 


0.127 


0.212 


1,491 


0.093 


0.05 


Ph-E-2 


1385.78 


0.65 


0.91 


3.48 


5.19 


0.067 


0.135 


0.149 


33,678 


0.101 


0.04 


Ph-E-4 


1399.96 


0.68 


0.91 


4.02 


4.82 


0.048 


0.079 


0.181 


1,547 


0.070 


0.04 


Ph-F-2 


1543.27 


0.60 


1.1 


3.81 


4.61 


0.053 


0.048 


0.186 


31,247 


0.095 


0.05 


Ph-F-4 


1559.44 


0.62 


1.1 


4.00 


4.54 


0.059 


0.057 


0.203 


1,547 


.075 


0.05 


Ph-G-2 


1561.75 


1.06 


0.18 


0.78 


3.33 


0.100 


0.221 


0.097 


42,528 


0.168 


0.17 


Ph-G-4 


1599.17 


1.04 


0.18 


1.10 


2.98 


0.109 


0.164 


0.116 


1,586 


0.140 


0.17 


Ph-H-2 


1676.43 


1.14 


0.21 


1.98 


4.66 


0.155 


0.212 


0.117 


35,048 


0.095 


0.1 


Ph-H-4 


1710.19 


1.14 


0.21 


2.75 


3.59 


0.109 


0.115 


0.178 


1,437 


0.069 


0.1 


Ph-I-2 


2236.05 


1.03 


0.56 


4.18 


4.86 


0.041 


0.059 


0.190 


35,754 


0.102 


0.02 


Ph-I-4 


2269.09 


1.05 


0.56 


4.48 


5.02 


0.045 


0.051 


0.208 


1,641 


0.073 


0.02 



Table 2. Basic structural parameters of Phoenix clusters at z = 0. The leftmost column labels each run, as in Table [T] the second and third colu mns list the 
peak circular v elocity, V ma x, and the radius, r max , at which it is reached. The concentration parameters of the best NFW tN avarro e t al. 1996, 1997) and Einasto 
lEinastoll 19651) fits are listed under c and ce, respectively. Qnfw and Qe are the figures of merit of the best NFW and Einasto fits, respectively. The column 
labelled a lists the Einasto shape parameter. N su ^ denotes the total number of subhaloes with more than 20 particles identified within T2oo; /sub is 'he fraction 
of the virial mass contributed by such subhaloes; and d B {{ is the distance from the gravitational potential minimum to the centre of mass of particles within the 
virial radius, in units of T2qq- 
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Name 


<Zh> 




</sub> 


(a) 


(Gmm) 


(Cx) 


(X) 


(Nm) 


<*> 




(d) 


Phoenix 


0.65 


0.061 


0.109 


0.175 


0.086 


1.75 


-1.86 


7866 


-0.97 


3984 


-3.32 




±0.36 


±0.047 


±0.027 


±0.046 


±0.041 


±0.29 


±0.04 


±965 


±0.02 


±317 


±0.10 


Aquarius 


1.65 


0.032 


0.071 


0.159 


0.048 


2.19 


-1.82 


5092 


-0.94 


4033 


-3.13 




±0.65 


±0.011 


±0.022 


±0.022 


±0.012 


±0.14 


±0.02 


±677 


±0.02 


±500 


±0.09 



Table 3. Comparison of the average properties of the six galaxy-sized Aquarius haloes and the nine cluster-sized Phoenix haloes. Sample averages are listed 
for each quantity together with the rms dispersion around the mean. The first column identifies the simulation set; zi, is the half-mass formation redshift; d a ff 
and / su (, are the dynamical relaxation diagnostics introduced in Table[2] a is the best-fit Einasto shape parameter and Q mm the goodness of fit measure (Sec. [6); 
Ci and x are the parameters of power-law fits to the pseudo-phase-space density profile (eq. [5J; N m and s describe the power-law fits to the subhalo mass 
function, N(> fj) — N m (fj/10~ 6 ) s (eq.[4); N v and d those corresponding to fits of the form, N(> v) = N r (v/0.025)'', to the subhalo velocity function (eq.|5J. 
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Figure 5. Spherically-averaged density (left panels) and logarithmic slope (right panels) of all level-2 Phoenix haloes as a function of radius. Radii have been 
scaled to the virial radius of each halo in the top panels and to the "scale radius", r_2, of the best-fit Einasto profile in the bottom panels. Profiles are plotted 
down to the convergence radius, r conv . The thick dashed black line shows the average density profile of all Phoenix haloes, computed after stacking the nine 
haloes, each scaled to its own virial mass and radius. The thick red dashed line shows the result of the same stacking procedure, but applied to the Aquarius 
haloes. 



regions where the two body-relaxation time exceeds the age of 
the Universe. This constraint defines a "convergence radius", r eonv , 
outside which the circular velocity, V c = (GM(< rf/r) 1 / 2 , is ex- 
pected to converge to better than 10%. Since V c is a cumulative 
measure we expect r conv to be a conservative indicator of the inner- 
most radius where local estimates of the density, p(r), converge to 
better than 10%. 

This is indeed the case for Ph-A, as shown in Fig. [4] The left 
panel shows p(r), multiplied by r 2 in order to remove the domi- 



nant radial trend so as to enhance the dynamic range of the plot. 
Thick lines highlight the radial range of the profile outside the con- 
vergence radius; the density clearly converges to better than 10% at 
radii greater than r com . In those regions the logarithmic slope yis 
also robustly and accurately determined. We conclude that r > r conv 
is a simple and useful prescription that identifies the regions unaf- 
fected by numerical limitations. We list r con v for all Phoenix runs 
in Table[[] 

The thin dashed lines in Fig|4] indicate the best-fit NFW 
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(brown) and Einasto (magenta) profiles, computed as described in 
the Appendix. The NFW shape is fixed in this log-log plot, whereas 
the Einasto shape is controlled by the parameter a, which is found 
to be 0.215 by the fitting procedure when applied to the Ph-A-1 pro- 
file. This figure suggests that the shape of the mass profile deviates 
slightly but systematically from the NFW profile. Although it is 
possible to obtain excellent fits over the resolved radial range with 
the NFW formula (typical residuals are less than ~ 10%), there is a 
clear indication that the density profile near the centre is shallower 
than the asymptotic r - 1 NFW cusp. In a greement with results from 
the Aquarius Project ( lNavarroetal][2 010). there is little indication 
that the central density cusp of Ph-A is approaching a power-law; 
the profile becomes gradually shallower all the way in to the inner- 
most resolved radius. This radial trend is very well described by the 
Einasto profile. 

Fig.[5]shows the density profiles of all level-2 Phoenix haloes, 
in a format similar to that of Fig. [4] The top panels show profiles 
with radii scaled to the virial radius of each cluster, whereas those 
at the bottom show radii in units of the "scale radius", r_2, of the 
best Einasto fit. Profiles are shown from the convergence radius, 
'conv, outwards. 

In general, the density profiles of Phoenix clusters become 
gradually shallower towards the centre; from y ~ 3 in the outer re- 
gions to an average value of y ~ 1 at the innermost resolved radius. 
This behaviour is similar to that of Aquarius haloes, whose average 
profile is shown by the thick red dashed lines in Fig. [5] The large 
difference between Aquarius and Phoenix seen in the top panels of 
this figure just reflect the different concentratior0 of cluster- and 
galaxy-sized haloes. Indeed, when radii are scaled to r_2, the aver- 
age Phoenix and Aquarius profiles are basically indistinguishable 
from each other. 

This is confirmed quantitatively by the best-fit Einasto pa- 
rameters of these average profiles (listed in Table [3]l. The average 
Phoenix halo is only slightly worse fit by an Einasto profile than 
Aquarius, as shown by the <2min goodness-of-fit measure (6.5% vs 
1.8%, respectively). There is also a slight difference in shape pa- 
rameter; the average Phoenix cluster has a = 0.175 whereas the 
average Aquariu s halo has a = . 159, in agreement with previously 
reported trends jGao et alj|2008l) . 

One aspect in which Phoenix and Aquarius haloes do differ is 
the halo-to-halo scatter: the dispersion in the Einasto parameter a 
is twice as large for clusters as for galaxy-sized haloes (Table [3](- 
This may be readily seen in Fig. [5] Ph-A-2, for example, follows 
the steady decline in y towards the centre exhibited by Ph-A- 1 (and 
characteristic also of Aquarius haloes), whereas in other cases, such 
as Ph-H-2, y stays roughly constant over a wide radial range near 
the centre. 

The latter behaviour is poorly captured by the Einasto or NFW 
fitting formulae, and leads to larger residuals and figure-of-merit 
values for the best fits. NFW and Einasto best-fit residuals are 
shown in Fig. |6j per bin deviations of up to 40% from NFW and 
~ 20% from Einasto fits are not uncommon for Phoenix clusters. 

These deviations may be traced to transient departures from 
equilibrium induced by the recent formation of many Phoenix clus- 
ters. For example, one of the worst offenders is Ph-H-2, which ac- 
creted half its final mass since z = 0.21 and whose unrelaxed ap- 
pearance is obvious in Fig. [3] In contrast, Ph-A-2, the cluster with 



2 The concentration is defined as r2oo / r~z, where r_2 is the radius at which 
the logarithmic slope y has the isothermal value of 2. This indicates the 



highest formation redshift of the Phoenix series (zj, = 1.17) is very 
well fit by both the Einasto and NFW profiles, with average resid- 
uals of only ~ 3% and ~ 6%, respectively. Indeed, a well defined 
correlation may be seen in Table [2] between quantitative measures 
of departures from equilibrium, such as the centre offset, d s, or 
the mass fraction in the form of substructure, / su (,, and the aver- 
age residuals from the best NFW and Einasto fits. On average, both 
indicators are substantially smaller for Aquarius than for Phoenix 
(Table [3}, as expected. The higher formation redshift of galaxy- 
sized haloes means that they are closer to dynamical equilibrium 
than recently-assembled cluster haloes. 



3.2 Pseudo-Phase-Space Density and Velocity Anisotropy 

The similarity in the mass profiles of galaxy- and cluster-sized 
CDM haloes highlighted in the previous subsection extends to their 
dynamical properties. We show this by comparing the spherically 
averaged pseudo-phase-space density (PPSD) profiles of Phoenix 
and Aquarius haloes. The PPSD, p/<7 3 , is dimensionally identi- 
cal to the phase-space density, but not strictly a measure of it. 
(Here the velocity dispersion, a(r), is defined as the square root 
of twice the specific kinetic energy in each spherical shell.) It is 
well known that PPSD profiles are well approx imated by a simple 
power- law, p/o 3 oc jTavlor & NavarrcfeoOlh . intriguingly sim- 



ilar to the secondary-infall self-similar solutions of Bertschinger 
(1985), where the exponent, x — —1-875 (see also iLudlow et al. 
2010, and references therein). 

PPSD profiles for all level-2 Phoenix clusters are shown in 
Fig. [7] and compared with the average PPSD for Aquarius haloes. 
Since clusters are denser and have higher velocity dispersions than 
galaxy-sized haloes, we scale all profiles to the scale radius, r_2, 
of each halo. Together with the density at this characteristic ra- 
dius, p_2. these quantities define a characteristic velocity, V_2 = 



(Gp_ 



U/2 



r_2, that allows us to compare PPSD profiles of haloes 



of widely different mass in a meaningful way. 

The top panel of Fig. UJ shows that, in these scaled units, 
Aquarius and Phoenix have very similar PPSD profiles. The simi- 
larity extends over the range 0.06 < r/r_2 < 4 where both simula- 
tion sets give converged results. (Note that Phoenix profiles actually 
probe radii interior to 0.06 r_2 because of their lower concentra- 
tion.) Table|3]lists the average parameters (and their dispersion) of 
power-law fits of the form 



P-2 



''-2 



(3) 



location of the maximum of the curves shown in Fig. [5] 



where C x = (c(r_2)/V_2) 3 - Fits are carried out over the length 
range 0.06 < r/r_2 < 4 for each halo. On average, both the slope 
and the normalisation of Aquarius and Phoenix PPSD profiles are 
almost indistinguishable emphasizing again the structural similar- 
ity between cluster- and galaxy-sized haloes. 

At the same time, the scatter is larger in the Phoenix sample 
than in Aquarius (Table highlighting again the larger halo-to- 
halo variation of cluster profiles. This may also be appreciated in 
the bottom panel of Fig. [7j where residuals from the self-similar 
r - 1.875 p 0wer i aw are shown. Although PPSD profiles scatter above 
and below the self-similar solution depending on the individual dy- 
namical state of each cluster, the PPSD profile of cluster haloes 
seem to be, on average, indistinguishable from that of galaxy-sized 
haloes. 

We reach a similar conclusion when comparing the veloc- 
ity anisotropy profiles of Phoenix clusters with those of Aquarius 
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Figure 6. Residuals from the best Einasto (left panel) and NFW (right panel) profile fits for all level-2 Phoenix haloes. Colors and line types are as in Fig[5] 
The thick black dashed curve corresponds to the composite profile obtained after stacking all 9 Phoenix level-2 runs. The red thick dashed curve corresponds 
to the same composite profile, but for the 6 galaxy-sized level-2 Aquarius haloes. 
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Figure 7. Spherically averaged pseudo-phase-space density (PPSD; p/o" 3 ) 
of all level-2 Phoenix haloes as a function of radius. Profiles are plotted 
down to the convergence radius, r conv . Radii are given in units of the scale 
radius, r_2> of the best-fit Einasto profile for each halo. Densities are scaled 
to p_2 = p(f-2), and velocity dispersions, a(r), to the characteristic veloc- 
ity V-2 = (Gp-2) 1 ^ 2 r-2- The thick dashed lines shows the average PPSD 
of all Phoenix (black) and Aquarius (red) haloes plotted over the converged 
radial range common to both simulation series: 0.06 < [r/r— 2) < 4, respec- 
tively. The bottom panel shows residuals from a simple r -1 ' 875 power-law 
fit. 



haloes (Fig. [8). Aside from a slightly larger scatter, the velocity 
anisotropy, which measures the ratio of the kinetic energy in tan- 
gential and radial motions, increases gently from the centre, where 
haloes are nearly isotropic, to the outer regions, where radial mo- 
tions dominate. Phoenix and Aquarius haloes again seem indistin- 
guishable from each other regarding velocity anisotropy when com- 
pared over their converged radial range. 



3.3 Projected Profiles 

The preceding discussion highlights the mass invariance of the 
structure of CDM haloes, but it also makes clear that the dynamical 
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Figure 8. Velocity anisotropy profiles of all level-2 Phoenix haloes. Radii 
are expressed in units of the scale radius, r-2, of the best-fit Einasto profile. 
Profiles are plotted down to the convergence radius, r conv . The thick dashed 
lines show the average anisotropy profile of all Phoenix (black) and Aquar- 
ius (red) haloes over the radial range where both give converged results 
(0.06 < {r/r - 2 ) < 4). 



youth of clusters limits the validity of simple fitting formulae to de- 
scribe their instantaneous mass profiles. This complication must be 
taken into account when comparing observational estimates of in- 
dividual cluster mass profiles with the profiles expected in a CDM- 
dominated Universe. Stacking clusters in order to obtain an "av- 
erage" cluster profile might offer a way of circumventing this dif- 
ficulty. This should smooth out local inhomogeneities in the mass 
distribution and average over different dynamical states to produce 
a more robust measure of the shape of the mass profile. 

Aside from dynamical youth, another issue that complicates 
the interpretation of observations is the fact that, due to the clus- 
ter's asphericity, projected mass profiles, such as those measured 
through gravitational lensing, may differ substantially from sim- 
ple projections of the 3D spherically-averaged profiles discussed 
above. 

Depending on the line of sight a cluster may appear more or 
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less massive within a given radius, leading to biases in the cluster's 
estimated mass, concentration, and even the shape of its density 
profile. We illustrate this in Fig. [9] where we plot the surface density 
profile of two Phoenix clusters, Ph-A-2 and Ph-I-2, each projected 
along 20 different random lines of sight. The aspherical nature of 
the clusters results in large variations (up to a factor of 3) in the 
surface density in the inner regions. For comparison, we also show 
in Fig.[9]the result of a wea k and strong-lensing analysis of a stack 
of four massive clusters bv lUmetsu et alj J201 lh . The mass of the 
stacked cluster lies between that of Ph-A and Ph-I, which explains 
why, on average, Ph-A T,(R) profiles lie below the observed data 
whereas the opposite applies to Ph-I. 

In agreement with earlier work (se e, e.g., ICorless & Kind 
l2007tlSereno et al.l201(jl ; IOguri et al.l201Ct and references therein), 
Fig. |9] suggests that substantial biases may be introduced by pro- 
jection effects on estimates of cluster parameters, especially when 
reliable data are restricted to the inner regions of a cluster. For 
example, fitting the inner 500/; kpc of the Ph-A-2 projected 
profile with an NFW profile results in mass and concentration 
(M200, c) estimates that vary from (5.4 x 10 14 A _1 Mq, 4.8) to 
(7.3 x 10 14 /i _1 M©, 9.8) when using the projections that maximize 
or minimize the inner surface density, respectively (see Fig. [9^. 
The corresponding numbers for Ph-I-2 are (1.8 x 10 15 h Mq, 4.1) 
and (3.0 x 10 15 h Mq, 7.1). Comparing these numbers with those 
listed in Tables[T]and[2]we see that variations as large as ~ 30% in 
the mass and ~ 60% in the concentration may be introduced just 
by projection effect^- 

We explore this further in Fig. [TO] where the small dots show 
the mass-concentration estimates for 500 random projections of 
each level-2 Phoenix cluster. Large symbols correspond to the 3D 
estimates listed in Tables[T]and[2] The black diamond symbol indi- 
cates the M?o n-c estimate for the s tack of 4 strong-lensing clusters 
presented by lUmetsu et alj J201 lh - This figure again emphasizes 
the importance of projection effects; for example, 12% of random 
projections result in concentration overestimates larger than 25%. 
Although an exhaustive analysis of such biases is beyond the scope 
of the present paper, the results in Figs.[9]and [l0]suggest that there 
is no substantial difficulty matching the surface densi t y pro file of 
lensing clusters such as those studied bv lUmetsu et al.1 d201 lh . Our 
interpretation thu s agrees with that reached by a number of recent 
studie s (see. e.gJOguri et alj201ll ; lOkabe et alfeoifj ; iGralla et all 



1201 it lUmetsu et al 2011 ), which conclude that there is no obvi- 
ous conflict between the concentration of lensing-selected clusters 
and those of ACDM haloes once projection effects are taken into 
account. Interestingly, despite the large variations in surface den- 
sity alluded to above, the shape of the surface density profile is 
quite insensitive to projection effects. We show this in the right- 
hand panels of Fig. [9] the weak dependence of y p (R) on projection 
may thus be profitably used to assess the consistency of theoretical 
predictions with cluster mass profiles. For illustration, we compare 
in the same panels the logarithmic slope of the proj ected profile, 
y g = d \nL(R) /dlnR, with the stacked cluster data of lUmetsu et all 
1 201 lh . Despite the fact that the mass of the simulated and observed 
clusters are different and that no scaling has been applied, there 
is clearly quite good agreement between observation and Phoenix 
clusters, supporting our earlier conclusion. Available data on indi- 
vidual clusters are bound to improve dramatically with the advent 



3 Note that variations may actually be larger, because these estimates ne- 
glect the possible contribution of the large-scale mass distribution along the 



Ph-A-2 Ph-E-2 
Ph-B-2 Ph-F-2 
Ph-C-2 Ph-G-2 Ph-I-2 
Ph-D-2 Ph-H-2 
J , , 1 I , 



♦ Umetsu et al. (2011) 
, I , , , L 



14.8 



15.0 

log10(M„h/M o ) 



15.2 



15.4 



Figure 10. Cluster virial mass vs concentration estimated from fits to the 
projected density profiles of level-2 Phoenix haloes in the radial range 
R < 500ft - 1 kpc. A total of 500 random projections are used for each halo. 
The large filled circles indicate the true value of the virial mass and concen- 
tration of the cluster, obtained from NFW fits to the 3D spherically-averaged 
profile (see Appendix and Table [2). The dashed curve flanke d by dotted 
lines s hows the fit to the mass-concentration relation derived bv lNeto et alj 
<2007l) . Note that projection effects lead to significant bias in the mass and 
concentration. , which are underestimated on average by 8.5 ± 17% and 
0.4 ±20%, respectively, where the "error" is the rms of all projections for 
the 9 clusters. The black diamond symbol indicat es the M200 -C estimate for 
a stack of 4 strong-lensing clusters taken from lUmetsu et al] i201lh . 



of surveys such as CLASH with the A dvanced Camera for Surveys 
onboard the Hubble Space Telescope jPostman et alj|20lll) . These 
surveys will enable better constraints on the shape of the inner mass 
profile of individual rich clusters and it is therefore important to 
constrain how projection effects may affect them. Fig. II llshows the 
distribution of y p at two projected radii, R = 3 and 10A kpc. The 
histograms are computed after choosing 500 random lines of sight 
for each of our 9 level-2 Phoenix haloes. On average, cluster pro- 
jected profiles flatten steadily towards the centre, from (y p ) = 0.35 
to 0.25 in that radial range, but with fairly large dispersion; the rms 
is Gy p = 0.054 and 0.091 at R = 10 and 3h~ [ kpc respectively. Be- 
cause of the large dispersion it is unlikely that observations of a 
single cluster can lead to conclusive statements about the viabil- 
ity of CDM; however, it should be possible to use this constraint 
fruitfully once data for a statistically significant number of clusters 
become available. 



line-of- sight. 



4 THE SUBSTRUCTURE OF PHOENIX CLUSTERS 

As may be seen from the images presented in Fig. [3] substruc- 
ture is ubiquitous in Phoenix clusters. We have used SUBFIND 
dSpringel et ail l2001bh to identify and characterize self-bound 
structures (subhaloes) within the virial radius of the main halo. 
We discuss below the mass function, spatial distribution, and in- 
ternal properties of subhaloes in Phoenix. Since our main goal is 
to explore the mass invariance of the properties of CDM haloes, 
we contrast these results with those obtained for the galaxy-sized 
Aquarius haloes. 
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Figure 9. Projected density profiles of Ph-A-2 (top) and Ph-I-2 (bottom). We show 20 different random projections for each cluster. The asphericity of the 
clusters leads to large variations (up to a factor of 3) in the projected density at a given radius depending on the line of sight. On the other hand, the shape 
of the profile (as measured by the logarithmic slope, y„ = — rflnZ/rflnR, is much l ess sensitive to pro jection effects. Data with error bars correspond to the 
stacked profile of 4 massive clusters estimated using strong and weak lensing data lUmetsu et al J201 ll) . 
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Figure 12. Left; The cumulative mass function of substructure haloes ("subhaloes") within the virial radius of cluster Ph-A at z = 0. We compare the results 
of four different realizations of the same halo, Ph-A-1 to Ph-A-4, with varying numerical resolution. The top and bottom panels contain the same information; 
the bottom shows the number of subhaloes weighted by mass or, equivalently, the fractional contribution of each logarithmic mass bin to the total mass in 
subhaloes. Each curve extends down to a mass corresponding to 60 particles. Note that, over the range resolved by the simulations, the cumulative function is 
well approximated by a power-law, N <* M _1 , the critical dependence for logarithmically divergent substructure mass. Right: Same as left panels, but for the 
subhalo peak circular velocity. 



The Phoenix project 13 



0.3 



0.2 



0.1 





I _ 

3 kpc/h I 






10 kpc/h J 










i , "i= — i 



0.0 0.2 0.4 0.6 



y p = -dlnL/dlnr 

Figure 11. Distribution of the slope of the circularly-averaged surface den- 
sity profile, y p (R), measured at two different radii, R = 3 and Wh kpc 
in projection. These histograms are based on 500 random lines of sight for 
each of the level-2 Phoenix clusters. Vertical arrows show the values corre- 
sponding to the projected profile of all nine clusters stacked together. The 
profiles become gradually shallower towards the centre, but with large scat- 
ter: (y p ) changes from 0.35 to 0.25 as R goes from 10 to 3ft kpc, but the 
halo-to-halo scatter is quite large, with rms of order 0.09 at 3 ft -1 kpc and 
0.05 at 10/i~' kpc, respectively. 
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Figure 14. Peak circular velocity, V max , vs the radius at which it is reached, 
'"max- The solid cyan curve indicates the r mM - y max relation obta ined for iso- 
lated haloes in the Millennium Simulation bv lNeto et alj j2007l) . Subhaloes 
in both Phoenix (solid black curve) and Aquarius (solid red curve) deviate 
systematically from this relation towards smaller r max at a given velocity. 
This is a result of tidal stripping, which shifts the location of the peak in- 
wards while changing little the peak velocity. Isolated haloes identified in 
Aquarius and Phoenix (shown with dashed lines) are not subject to tides 
and are in good agreement with the Millennium Simulation results. 



4.1 Mass Function 

We start by analyzing the Ph-A simulation series in order to iden- 
tify the limitations introduced by finite numerical resolution. The 
top left panel of Fig. [12] shows the cumulative mass function of 
subhaloes, N(> M), plotted in each case down to the mass cor- 
responding to 60 particles. The bottom left panel shows the same 
data, but after weighting the numbers by subhalo mass, M su b, in 
order to emphasize the differences between runs. The results show 
clearly how, as resolution improves, the mass function converges at 
the low-mass end. Ph-A-4 agrees with higher resolution runs for 
subhaloes with mass exceeding ~ 2 x W Mq, correspond- 
ing to roughly 150 particles; the same applies to Ph-A-3 for mass 
greater than ~3x 10 9 hT x M Q , or ~ 170 particles, and to Ph-A-2 
for ~7x 10 8 h~ l M , or 140 particles. We conclude that the sub- 
halo mass function can be robustly determined in Phoenix haloes 
down to subhaloes containing roughly 150 particles, in good agree- 
ment with the result s reported for Aquarius haloes (see Fig. 6 of 
ISpringel et ail20 08a). For level-2 runs this implies a subhalo mass 
function that spans over 6 decades in mass below the virial mass 
of the halo. The subhalo mass function is also routinely expressed 
in terms of the subhalo peak circular velocity. This is shown in 
the right-hand panels of Fig. [T2] which shows that level-2 Phoenix 
runs give robust estimates of the abundance of subhaloes down to 
Vmax ~ 20 km s~ 1 , a factor of ~ 75 lower than the main halo's V200 • 
Both the subhalo mass and velocity functions seem reasonably 
well approximated by simple power laws: N and N <x V^^ 4 , 

respectively. Interestingly, the M~ 1 dependence corresponds to the 
critical case where each logarithmic mass bin contributes equally 
to the total mass in substructure. This is logarithmically divergent 
as M su [, approaches zero, and implies that a significant fraction of 
the mass could in principle be locked in haloes too small to be 
resolved by our simulations. We note, however, that even at the 
resolution of Ph-A-1, of nearly 7 decades in mass, only 8% of the 
mass within r2oo is in the form of substructure. Extrapolating down 



to the Earth mass by assuming that N °c M~ ub , the total mass locked 
in substructure would still be only about 27 percent. 

Fig.Q~3]compares these results with other level-2 Phoenix clus- 
ters in order to assess the general applicability of the Ph-A subhalo 
mass function. The cumulative number of subhaloes N(> M) is 
weighted here by /j = M su b/M2oo (left panel) in order to emphasize 
differences as well as to enable the comparison of haloes of dif- 
ferent virial mass. Although the subhalo mass function, expressed 
in this form, is relatively flat in several Phoenix clusters (indica- 
tive of an N « A^sub dependence) it is clearly declining in others. 
The average trend, as indicated by the "stacked" Phoenix clus- 
ter (thick dashed black curve) may be approximated, in the range 
10~ 6 < u < 10~ 4 , by N <x ^ _0 - 98 . This is a slightly steeper depen- 
dence than found for Aquarius haloes over the same mass range, 
N 0= jU -0,94 (thick dashed red curve), but still subcritical. The slight 
difference in the average slope of the Aquarius and Phoenix sub- 
halo mass functions is smaller than the halo-to-halo scatter in either 
simulation set. This is shown in Table [3] where we list the average 
parameters of power-law fits of the form, 

N(>u)=N m (p/10- 6 ) s , (4) 

for Aquarius and Phoenix haloes. The dispersion around (s) is sim- 
ilar to the difference between the average slope of Aquarius and 
Phoenix haloes, suggesting that there is no significant difference in 
the shape of the subhalo mass function of cluster- and Milky-way 
halo-sized haloes. 

Fig. [T3] also shows that substructure is slightly more preva- 
lent in clusters than in galaxy-sized haloes. Indeed, at all values 
°f ^sub/^200 me number of Phoenix subhaloes exceeds that in 
Aquarius, and this is reflected in the higher values of (N m ) (7866 
for Phoenix vs 5092 for Aquarius; see TablefJJ. This is another con- 
sequence of the dynamical youth of clusters compared to galaxies 
(tides take a few orbital times to strip a subhalo), as may be verified 
by inspection of TableQ] in the cluster that forms latest, Ph-G, sub- 
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Figure 13. As the bottom panels of Fig. 1121 but for all level-2 Phoenix haloes. The cumulative mass function (left panel) is weighted by subhalo mass, expressed 
in units of the virial mass. A cumulative N « M~' dependence, the critical case for logarithmic divergence in the total substructure mass, corresponds to a 
horizontal curve in these scaled units. Although the dependence is nearly flat in several Phoenix clusters it is clearly declining in others, and the average trend 
seems to be sub-critical. Compared with Aquarius (thick dashed red curve) the average Phoenix subhalo mass function is slightly steeper. The panel on the 
right is analogous to the mass function, but for the subhalo peak circular velocity, weighted by V,^ ax . (See text for further discussion.) 



structure makes up roughly 17% of its virial mass, but in the earliest 
collapsing system of the Phoenix series, Ph-A, it makes only 8%. 

Interestingly, as a function of v = Vm ax /V200, the compari- 
son between the Aquarius and Phoenix subhalo functions reverses 
(right-hand panel of Fig. 1 1 3b - At a given velocity (scaled to the 
virial value), subhaloes are more abundant in Aquarius than in 
Phoenix. This is a consequence of tidal stripping, which affects 
Aquarius subhaloes more: since tides act to remove preferentially 
the outer regions of a subhalo, they affect more its mass than its 
peak circular velocity. 

For example, as discussed by IPenarrubia et alj ([2008), after 
losing half of its mass to tides, the peak velocity of a subhalo de- 
creases only by ~ 25%. Even after losing 90% of its mass, V max is 
only reduced by about one half. Since Aquarius haloes form earlier, 
their subhaloes have been accreted earlier and, on average, have 
been more stripped than Phoenix subhalos, leading to higher rela- 
tive velocities for their bound mass than for Phoenix subhalos. This 
shifts their abundance when measured in terms of peak velocity. In 
the range 0.025 < V < 0.1 fits to the subhalo function of the form 

N{>\) = M, (v/0.025) d (5) 

yield (N v ) = 4033 and (d) = -3.13 for Aquarius and 3984 and 
—3.32, respectively, for Phoenix (see Table [3). Given the scatter, 
the difference seems too small to be significant. We conclude that 
the scaled subhal o velocity functio n, N(>v), is roughly indepen- 
dent of mass (see lWang et alj2012l. for a more thorough discussion 
of this point). 

The effects of tidal stripping on Phoenix subhaloes is shown 
in Fig. 1 141 Here we plot Vmax vs r max for subhaloes identified in 
Ph-A-1 (solid black curve). This relation is clearly offset from the 
mean relation that h olds for iso l ated h aloes in the Millennium Sim- 
ulation, as given by |NetoetalH2007h (cyan line). As expected for 
haloes that have undergone tidal stripping, r max shifts inwards as 
the subhalo loses mass whilst leav ing the peak velocity relatively 
unchanged I Pen arrubia et alj[2008h . Support for this interpretation 
may be found by inspecting the same relation for "isolated" haloes 
in Phoenix (i.e., those outside the main halo and that are not em- 
bedded in a more massive structure; the r max -V max relation for these 



systems (see dashed lines) is consistent with that of Millennium 
haloes. 

Fig. [11] also includes results for isolated haloes and subhaloes 
in Aquarius (red lines). The results from the two sets of simula- 
tions form a single sequence and this allows us to characterize the 
structural parameters of subhaloes over a range spanning more than 
two decades in velocity (and thus over six decades in mass). On 
average, subhaloes follow the same r max -V max scaling relations as 
isolated haloes, but shifted by about a factor of two in radius (or, 
alternatively, by ~ 30% in velocity). 

We conclude from this discussion that although substructure 
does not seem fully invariant with halo mass, the changes are rel- 
atively small when comparing the haloes of clusters and galaxies, 
and depend on whether subhalo masses or velocities are used to 
characterize substructure. The subhalo mass function of clusters, 
scaled to halo virial mass, is similar in shape to that of galaxy-sized 
haloes (which are roughly one thousand times less massive), but 
with a slightly higher normalization (~ 35%). The normalization 
difference disappears when the scaled subhalo velocity function, 
N(> v), is used. The total mass in substructure increases with the 
dynamical youth of the system and is more prevalent in clusters 
than on galaxy scales, but only weakly so: the average mass frac- 
tion in substructures is 1 1% for Phoenix and 7% for Aquarius. 



4.2 Spatial Distribution 

The distribution of subhaloes within the main ha l o has been the 
subjec t of many studies (e.g . iGhigna et al.l |2000J: iDiemand et al. 
2004b:; |Pe Lucia et alj 12004 ; bap et alj l2004bUiJ; ISpringel et all 
2008a; Ludlow et al. 2009) over the past decade. This work has 
demonstrated that substructure does not follow the same spatial 
distribution as the dark matter: subhaloes tend to populate prefer- 
entially the outskirts of the main halo and their spatial distribution 
is much more extended than the mass. It also hinted that the num- 
ber density profile of subhaloes is roughly independent of subhalo 
mass, at least in the subhalo mass range where simulations resolve 
them well and where they exist in sufficient numbers for their spa- 
tial distribution to be determined. This result has been confirmed 
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Figure 15. Subhalo number density profiles. The panel on the left shows the spatial distribution of subhaloes with more than 100 particles in each of the 9 
Phoenix level-2 clusters. Each profile is normalized to the mean number density of subhaloes within the virial radius. The thick dashed black curve traces 
the result of stacking all 9 level-2 Phoenix haloes. The profile obtained after stacking all level-2 Aquarius haloes is shown by the red dashed curve. Note that 
subhaloes are slightly more concentrated in the case of Phoenix than of Aquarius. The panel on the right shows the density profile of subhaloes in different 
bins of subhalo mass, computed after stacking all 9 level-2 Phoenix clusters. Note that the spatial distribution of subhaloes is approximately independent of 
subhalo mass. 



recently by the Aquarius simulati on suite for haloes similar to the 
Milky Way dSpringel et a"i]|2008ij) . 

A number of observational diagnostics depend on the spatial 
distribution of substructure, and it is therefore important to ver- 
ify that this result holds also on galaxy cluster scales. For exam- 
ple, recent analyses indicate that total flux of dark matter annihila- 
tion radiation is expected to be dominated by low-mass subh aloes 
dKuhlen et al1l2008l ; ISpringel et alj[2008bl : iGao et alj|201 lbh . It is 
therefore crucial to constrain their spatial distribution in order to 
understand the expected angular distribution of the annihilation 
flux and to design optimal filters t o aid its discovery (see, e.g., 
IPinzke et alj201 lhlGao et al.11201 lbh . 

We show the number density profile of subhaloes in Fig. [15] 
The left panel shows the profiles for each of the 9 level-2 Phoenix 
haloes (thin lines), as well as the profile corresponding to stacking 
all 9 haloes after scaling them to the virial mass and radius of each 
cluster (thick dashed black curve). All subhaloes with more than 
100 particles have been used for this plot. This figure clearly con- 
firms the results of earlier work: the subhalo distribution is more 
extended than that of the dark matter; In addition there is a well 
defined "core" in the central density of the subhalo distribution; 
Subhaloes primarily populate the outskirts of the main halo. 

There is also considerable halo-to-halo scatter, especially near 
the centre, where the number density of subhaloes may vary by up 
to a factor of three. Comparing the average number density profile 
of Phoenix with that of Aquarius (thick red dashed curve) reveals 
that cluster subhaloes are slightly more abundant near the centre, 
by up to 50% at r = 0.1r2oo- In me outskirts of the main halo 
both Aquarius and Phoenix give similar results. As discussed by 
iLudlow et al.l d2009r) , the number density profile can be fitted accu- 
rately by an Einasto profile (eq.|7J, just like the dark matter, but with 
quite different shape parameters: a ~ 1 for subhaloes but ~ 0.2 for 
the main halo. An Einasto fit to the Phoenix subhalo profile yields 
r_2 = 0.58 /^oo and a — 1-0. For Aquarius, the same procedure 
yields r_2 = 0.64r2oo and a = 1.0, and a central density normal- 
ization lower by a factor of 1.3, when expressed in units of («), the 
mean number density of subhaloes within r2oo- 
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Figure 16. Stacked subhalo number density profiles as a function of rj r2oo 
for the nine Phoenix haloes and for different definitions of the lower subhalo 
"mass" limit. The solid line shows the radial profile for all subhaloes whose 
progenitors had a maximum circular V max exceeding 45 km s _I when they 
first fell into the cluster; the dot-dashed line shows a similar profile but 
for subhaloes with V max greater than 30km s _I at the present day; finally, 
the dashed line show the profile for all subhaloes containing more than 200 
bound particles. For comparison, a dotted line shows the stacked dark matter 
mass profile of the clusters. The profiles are normalised to integrate to the 
same value within r^oo- Note that none of the subhalo profiles matches the 
shape of the dark matter profile within 0.25r2oo- 



with galaxies make a variety of assumptions about how to assign 
galaxies to subhaloes. A number of authors have argued that al- 
though present subhalo mass and maximum circular velocity are 
strongly affected by tidal stripping and so are poor indicators of 
galaxy properties, the mass or circular velocity at infall are plausi- 
bly much better and give mea ningful results when used in subhalo 
abundance matching analy s es (fVa le & Ostriker 2004; Conro v et al.l 
l2006l : lBehroozi et al.l20"Tfj|;lGuo et alj2010h . We study this issue in 
Fig. UU which shows stacked number density profiles for subhalo 
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samples defined above thresholds in present mass, present circular 
velocity and infall circular velocity. Note that these thresholds are 
chosen so that each sample contains roughly 6000 subhaloes. In 
agreement with earlier work, we see that sample definition has a 
substantial effect on the inferred radial profile of the subhalo pop- 
ulation. Subhalo samples defined by present mass have shallower 
profiles than samples defined by present circular velocity which, in 
turn, have shallower profiles than samples defined by infall circular 
velocity. Note, however, that all these profiles differ substantially 
from the mean dark matter density profile, especially in the inner 
regions (r < 0.25 rjoo), whereas observations show the mean galaxy 
number density profiles in the inner regions of clu sters to follow 
the mean dark matter profiles quite closely (e.g . ICarlberg et al.l 
ll997l ; lBiviano & Girardil2003l ; ISheldon et alj|2009l) . Semi-analytic 
models which explicitly follow the formation of galaxies within 
the evolving subhalo population provide a better match to the ob- 
served inner profiles because they include a population of "orphan" 
galaxies whose dark matter subhaloes have already been tida lly de- 
stroyed dGao et al j2004al : IWang e7ail2006l : lGuo et al.ll201 lb . 

Fig. [T7] shows the fractional contribution of substructure to 
the total mass of the halo, as a function of radius, either in cumula- 
tive (left panel) or differential (right panel) form. This figure shows 
quantitatively that substructure contributes only a small fraction of 
the halo mass. This contribution peaks in the outer regions; it is only 
0.1% at r = 0.02 r2oo but it reaches 10-20% at the virial radius. The 
total mass contribution is, on average, just over 10% (see also Ta- 
ble^- Results for Phoenix are similar to Aquarius, adjusted up by a 
modest amount that reflects the overall larger substructure fraction 
present in clusters relative to galaxy-sized haloes. This adjustment 
is mainly noticeable in the inner regions, reflecting our earlier con- 
clusion that substructure in Phoenix is more centrally concentrated 
than in Aquarius. 



5 SUMMARY AND CONCLUSIONS 

We present the Phoenix Project, a series of simulations of the 
formation of rich galaxy cluster halos in the ACDM cosmogony. 
Phoenix simulations follow the dark matter component of 9 dif- 
ferent clusters with numerical resolution comparable to that of 
the Milky Way-sized haloes targeted in the Aquarius Project 
dSpringel et"aT]|2008al : lNavarro et alj|2010h . We report here on the 
basic structural properties of the simulated clusters and their sub- 
structure, and compare them with those of Aquarius haloes in order 
to highlight the near mass invariance of cold dark matter haloes in 
the absence of baryonic effects. Our main results may be summa- 
rized as follows. 

Radial Profiles. The recent formation of galaxy clus- 
ters causes many of them to be rapidly evolving and unrelaxed. This 
results in mass profiles that are less well approximated by simple 
fitting formulae such as the NFW or Einasto profiles than those 
of galaxy haloes. Stacking clusters helps to average out inhomo- 
geneities in the mass distribution characteristic of transient states. 
The mass profile of the stacked cluster is very similar to that of 
Aquarius haloes; it can be well approximated by an Einasto pro- 
file, albeit with a slightly larger value of the shape parameter, a, 
and significantly lower concentration. The similarity extends to the 
dynamical properties of the haloes: when properly scaled, the aver- 
age pseudo-phase-space density and velocity anisotropy profiles of 
Aquarius and Phoenix haloes are indistinguishable. 

Density Cusp . The central density cusp has, at the inner- 
most resolved radius (r conv ~ 2 x 10~ 3 >"2ooX an average logarith- 



mic slope (y) = 1.05 ±0.19, where the "error" refers to the halo- 
to-halo rms dispersion of the 9 level-2 Phoenix runs. This is only 
slightly steeper than that of Aquarius haloes at comparable radii, 
for which (y) = 1.01 ±0.10). Although in some clusters y remains 
roughly constant over a sizeable radial range near the centre, in the 
majority of cases the profile keeps getting shallower all the way to 
the innermost converged radius, with little evidence of convergence 
to an asymptotic power-law behaviour. 

Projected Profiles. Because of their aspherical na- 
ture, the surface density of Phoenix haloes varies greatly depend- 
ing on the line of sight, in some cases by up to a factor of ~ 3 at 
a given projected radius. This affects especially the inner regions 
and may give rise to substantially biased estimates of a cluster's 
total mass and concentration. For example, NFW fits to the inner 
500 kpc of 9 Phoenix haloes, on average, lead to estimates of 
M200 and c that can be overestimated by 20% and 80%, respec- 
tively, when the cluster is projected along the major axis and under- 
estimated by 30% and 20% respectively when seen along the minor 
axis. The shape of the surface density profile, on the other hand, is 
hardly affected by projection. The average logarithmic slope of the 
surface density profile declines gradually towards the centre, from 
(j p ) =0.35 ±0.091 MR= 10/j" 1 kpc to 0.21 ±0.054 &tR = 3/T 1 
kpc, again with no clear sign of approaching a power-law asymp- 
totic behaviour. 

Substructure Mass Function . Substructure is more 
abundant (by about ~ 35% on average) in Phoenix clusters than in 
galaxy haloes. At a given /j = Af su b/Af2oo> the cumulative number 
of cluster subhaloes is higher in Phoenix by about ~ 30% com- 
pared to Aquarius, with a tendency for the excess to increase at the 
low-mass end. In some cases the subhalo mass function is best ap- 
proximated by a power law with the critical slope N °= . There 
is significant halo-to-halo scatter, however, and the average trend 
is subcritical. In the range 1 x 10~ 6 </)<lx 10~ 4 we find that 
N = 0.010/j -0,98 fits well the composite subhalo mass function of 
the 9 level-2 Phoenix clusters stacked together. For comparison, the 
same procedure for the Aquarius haloes yields a very similar result: 
Af = 0.012 J u" ' 94 . 

Substructure Spatial Distribution. We con- 
firm earlier reports that subhaloes are biased tracers of the halo 
mass distribution, avoiding the central regions and increasing in 
prevalence gradually from the centre outwards. As in galaxy haloes, 
the subhalo number density profile appears to be independent of 
subhalo mass, and may be approximated accurately by an Einasto 
profile, but with scale radius ~ 0.58 ^200 and a shape parameter 
much greater than that of the dark matter distribution, a ~ 1.0. 
Phoenix subhaloes are slightly more concentrated than those of 
Aquarius haloes: inside 0.1 ^200 they make up roughly 0.05% of 
the enclosed mass, a factor of 2 to 3 times larger than in Aquar- 
ius haloes. The difference decreases with increasing radius; in total 
Phoenix subhaloes make up on average 1 1% of the total mass, com- 
pared with 7% for Aquarius. 

Our analysis confirms the remarkable structural similarity of 
CDM haloes of different mass, whilst at the same time emphasiz- 
ing the small but systematic differences that arise as halo mass in- 
creases from galaxies to clusters. Many of these differences may 
be ascribed to the dynamical youth of galaxy clusters, which lead 
to larger deviations of individual clusters from the average trends. 
This argues for combining the results of as many clusters as pos- 
sible in order to average over the transient features of individual 
systems and to uncover robust trends that may be fruitfully com- 
pared with the predictions of the ACDM paradigm. 
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Figure 17. Left panel: Cumulative fractional contribution of subhaloes (resolved with more than 100 particles) to the enclosed mass, shown as a function 
of radius for all level-2 Phoenix clusters (thin lines). A thick dashed black curve shows the average trend, computed after stacking all 9 Phoenix haloes. 
The corresponding result for Aquarius is shown by the thick dashed red curve. Right panel: Fraction of total mass contributed by substructure in different 
radial bins. As in the left panel, only subhaloes with more than 100 particles are considered; black and red thick dashed lines correspond to the average trend 
computed after stacking all level-2 Phoenix and Aquarius haloes, respectively. 
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Because these formulae define the characteristic parameters in 
a slightly different way, we choose to reparametrise them in terms 
of r_2 and p_2 = p(r_2), which identify the "peak" of the r 2 p 
profile shown in the left panel of Fig.|4] This marks the radius where 
the logarithmic slope of the profile, y(r) = — dlnp/d\nr, equals 
the isothermal value, y=2. We note that, unlike NFW, when a is 
allowed to vary freely the Einasto profile is a 3-parameter fitting 
formula. 

6.2 Fitting procedure 

We compute the density profiles of each halo in 32 radial bins 
equally spaced in log 10 r, in the range r conv < r < r2oo- All haloes 
are centred at the minimum of the gravitational potential. Best-fit 
parameters are found by minimizing the deviation between model 
and simulation across all bins in a specified radial range. In the 
case of the density profile, the best fit is found by minimizing the 
figure-of-merit function, Q 2 , defined by 

2 2 =t— £(lnp,-ln P r de1 ) 2 . (8) 

■"bins ; = i 

This function provides a simple measure of the level of dis- 
agreement between simulated and model profiles. It is dimension- 
less; it weights different radii logarithmically; and, for a given ra- 
dial range, Q 2 is roughly independent of the number of bins used 
in the profile. The actual value of Q is thus a reliable and objective 
measure of the average per-bin deviation from a particular model. 
Thus, minimizing Q 2 yields for each halo well-defined estimates of 
a model's best-fit parameters. The values of 2min f° r eacn halo are 
given in Table|2]for both Einasto and NFW fits. 

It is less clear how to define a goodness-of-fit measure as- 
sociated with Q 2 and, consequently, how to assign statistically- 
meaningful confidence interva ls to the best-fit para meter values. 
We have explored this issue in iNavarro et ail £2010) and we refer 
the interested reader to that paper for details. 



6 APPENDIX 

6.1 Fitting formulae 

The fitting formulae used to describe the mas s profile of Phoenix 
haloe s are the following: (i) The NFW profile dNavarro et al.l 199(3 . 
1 19971) . given by 

{r/r s ){l + r/r s ) 1 
and (ii) the Einasto profile teinastdll965l) , 

ln(p(r)/p_2) = (-2/a)[(rA_2) a -l]. (7) 



